R Machine Learning - Kapitel 3
Einleitung zu Machine Learning
Machine Learning (ML) in der Statistik bezieht sich auf die Verwendung von statistischen Methoden und Algorithmen, um Modelle zu entwickeln, die es Computern ermöglichen, Muster in Daten zu erkennen und Vorhersagen zu treffen. Das Ziel von ML-Methoden besteht darin, aus Beispielen zu lernen, um Vorhersagen auf neuen, unbekannten Daten treffen zu können. Dies geschieht durch die Verwendung von Algorithmen wie Entscheidungsbaum-Modellen, künstlichen neuronalen Netzen, Bayes’schen Netzen und Support Vector Machines.
In der Praxis wird Machine Learning in der Statistik für viele Anwendungen eingesetzt, wie z.B. Bild- und Spracherkennung, Empfehlungssysteme, Betrugserkennung, medizinische Diagnose und Vorhersage von Aktienkursen. Die Wahl des geeigneten Algorithmus hängt von der Art der Daten, der Komplexität des Problems und den verfügbaren Ressourcen ab.
Supervised und unsupervised Modelle sind zwei grundlegende Kategorien von Machine-Learning-Modellen.
Supervised Learning
Supervised Learning ist eine Methode des maschinellen Lernens, bei der der Algorithmus trainiert wird, um aus einer gegebenen Menge von Trainingsdaten zu lernen, die annotiert sind, d.h. bei denen die Ausgabe oder Zielvariable bekannt ist. Der Algorithmus lernt, indem er die Beziehung zwischen Eingabe- und Ausgabevariablen modelliert, um Vorhersagen auf neuen Daten zu treffen. Es werden also Attribute X benutzt, um eine Zielvariable Y vorherzusagen. Die meisten Supervised Learning Probleme können einer von zwei Kategorien zugeordnet werden: Regression oder Klassifikation. Regressionsprobleme haben eine metrische Zielvariable, Klassifikationsprobleme haben eine kategoriale Variable.
Unsupervised Learning
Unsupervised Learning ist eine Methode des maschinellen Lernens, bei der der Algorithmus nicht mit annotierten Trainingsdaten arbeitet. Stattdessen lernt der Algorithmus, Muster und Strukturen in den Daten zu erkennen, indem er versucht, sie automatisch zu gruppieren oder zu segmentieren, um inhärente Beziehungen zu entdecken. Ein Beispiel für die Anwendung von Unsupervised Learning ist die Segmentierung von Kunden in Gruppen mit ähnlichem Kaufverhalten. Dabei werden keine annotierten Daten verwendet, sondern der Algorithmus erkennt die Gemeinsamkeiten in den Kaufmustern der Kunden und gruppiert sie entsprechend.
Bekannt sind hier vor allem das Clustering und die Dimensionsreduktion. Bei diesen Modellen ist die Qualitätsüberprüfung natürlich sehr schwierig, da wir Antwort auf unsere Frage nicht kennen. Im Gegensatz zu supervised Modellen.
Modellierungsverfahren
Sich den Machine Learning Modellen anzunähern bedeutet, die Daten besser kennezulernen. Validierungsprozess durchführen ist ebenso wichtig, wie die richtigen Features und Outcome Variablen zu wählen. So ist der Modellierungsprozess kein Sprint, sondern es bedarf vieler Schritte, um das optimale Modell zu finden.
Voraussetzungen
# Helper packages
library(dplyr) # for data manipulation## Warning: Paket 'dplyr' wurde unter R Version 4.2.3 erstellt##
## Attache Paket: 'dplyr'## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2) # for awesome graphics## Warning: Paket 'ggplot2' wurde unter R Version 4.2.3 erstellt# Modeling process packages
library(rsample) # for resampling procedures## Warning: Paket 'rsample' wurde unter R Version 4.2.3 erstelltlibrary(caret) # for resampling and model training## Warning: Paket 'caret' wurde unter R Version 4.2.3 erstellt## Lade nötiges Paket: latticelibrary(h2o) # for resampling and model training## Warning: Paket 'h2o' wurde unter R Version 4.2.3 erstellt##
## ----------------------------------------------------------------------
##
## Your next step is to start H2O:
## > h2o.init()
##
## For H2O package documentation, ask for help:
## > ??h2o
##
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit https://docs.h2o.ai
##
## ----------------------------------------------------------------------##
## Attache Paket: 'h2o'## Die folgenden Objekte sind maskiert von 'package:stats':
##
## cor, sd, var## Die folgenden Objekte sind maskiert von 'package:base':
##
## %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
## colnames<-, ifelse, is.character, is.factor, is.numeric, log,
## log10, log1p, log2, round, signif, trunc# h2o set-up
h2o.no_progress() # turn off h2o progress bars
h2o.init() # launch h2o## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 1 days 5 hours
## H2O cluster timezone: Europe/Berlin
## H2O data parsing timezone: UTC
## H2O cluster version: 3.40.0.1
## H2O cluster version age: 8 months and 1 day
## H2O cluster name: H2O_started_from_R_nikla_ijy204
## H2O cluster total nodes: 1
## H2O cluster total memory: 1.94 GB
## H2O cluster total cores: 4
## H2O cluster allowed cores: 4
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## R Version: R version 4.2.2 (2022-10-31 ucrt)## Warning in h2o.clusterInfo():
## Your H2O cluster version is (8 months and 1 day) old. There may be a newer version available.
## Please download and install the latest version from: https://h2o-release.s3.amazonaws.com/h2o/latest_stable.htmlWir arbeiten mit den Datensätzen Ames housing data und
Job attrition data.
library(modeldata)## Warning: Paket 'modeldata' wurde unter R Version 4.2.3 erstelltdim(modeldata::attrition)## [1] 1470 31churn <- attrition |>
mutate_if(is.ordered, .funs = factor, ordered = FALSE)
# Ames housing data
ames <- AmesHousing::make_ames()
ames.h2o <- as.h2o(ames)
# Job attrition data
#churn <- mutate_if(attrition, is.ordered, .funs = #factor, ordered = FALSE)
#churn.h2o <- as.h2o(churn)Daten splitten
Als Anforderung für unseren Algorithmus f(X) stellen wir, dass er akurat zukünftige Werte, die auf einer Menge von Features basieren, vorhersagt. Das nennt man Generalisierbarkeit unseres Algorithmus. Dafür müssen wir unseren Datensatz in zweit Teile teilen:
Trainingsdatensatz: Er entwickelt unsere Feature Menge, trainiert unseren Algorithmus, passt Hyperparameter an, vergleicht Modelle und entscheidet sich für ein finales.
Testdatensatz: Ist das finale Modell gewählt, so benutzen wir den Testdatensatz, um die Performance zu bewerten, die wir als Generalisierungsfehler bezeichnen.
Wir splitten die Trainingsdaten - Testdaten für gewöhnlich im Verhältnis 60 - 40, 70 - 30 oder 80 - 20 auf. Ist das Verhältnis zu groß, so finden wir ein Modell, dass die Trainingsdaten sehr gut anpasst, aber nicht generalisierbar ist (Overfitting). Die Testdaten sind möglicherweise nicht repräsentativ für die Gesamtpopulation. Ist der Testdatensatz zu groß, so können wertvolle Daten möglicherweise nicht in das Modell aufgenommen werden, da sie dem Testdatensatz zugeordnet sind und nicht im Trainingsdatensatz enthalten sind. Zu große Trainingsdatensätze führen auch irgendwann nicht mehr zu zustzlichem Informationsgewinn. Kleinere Sätze erhöhen auch die Geschwindigkeit.
Die zwei einfachsten Wege Daten in zwei zu teilen, sind die einfache Zufallsauswahl und das geschichtete Sampling.
Einfaches zufälliges Auswählen
# Using base R
set.seed(123) # for reproducibility
index_1 <- sample(1:nrow(ames), round(nrow(ames) * 0.7))
train_1 <- ames[index_1, ]
test_1 <- ames[-index_1, ]
# Using caret package
set.seed(123) # for reproducibility
index_2 <- createDataPartition(ames$Sale_Price, p = 0.7,
list = FALSE)
train_2 <- ames[index_2, ]
test_2 <- ames[-index_2, ]
# Using rsample package
set.seed(123) # for reproducibility
split_1 <- initial_split(ames, prop = 0.7)
train_3 <- training(split_1)
test_3 <- testing(split_1)
# Using h2o package
split_2 <- h2o.splitFrame(ames.h2o, ratios = 0.7,
seed = 123)
train_4 <- split_2[[1]]
test_4 <- split_2[[2]]Geschichtetes Sampling
Wir können aber auch einen Datensatz so auswählen, dass unsere Trainings- und Testdatensätze eine ähnliche Y Verteilung haben. Das bietet sich an, wenn wir z.B. bei Klassifizierungsproblemen eine sehr unbalanzierte Response Variable haben (in etwa 85% - 15%). Das rsample Packet bietet hier eine sehr unkomplizierte Lösung. Durch das Erzwingen von geschichtetem Sampling erzwingen wir hier Datensätze, die eine ähnliche Verteilung der Response haben.
# orginal response distribution
table(churn$Attrition) %>% prop.table()##
## No Yes
## 0.8387755 0.1612245##
## No Yes
## 0.8387755 0.1612245
# stratified sampling with the rsample package
set.seed(123)
split_strat <- initial_split(churn, prop = 0.7,
strata = "Attrition")
train_strat <- training(split_strat)
test_strat <- testing(split_strat)
# consistent response ratio between train & test
table(train_strat$Attrition) %>% prop.table()##
## No Yes
## 0.8394942 0.1605058##
## No Yes
## 0.838835 0.161165
table(test_strat$Attrition) %>% prop.table()##
## No Yes
## 0.8371041 0.1628959##
## No Yes
## 0.8386364 0.1613636Modelle erstellen
In R gibt es eine Vielzahl an leistungsstarken Algorithmen in verschiedensten Paketen. Manche sind rechnerisch effizienter als andere, manche sind flexibler. In den folgenden Kapiteln lernen wir verschieden Algorithmen und Pakete für unterschiedliche Problemstellungen kennen.
Resampling Methoden
Wir haben sehr explizit darauf hingewiesen, dass wir das Testset während der Trainingsphase nicht verwenden, um die Leistung des Modells zu bewerten. Wie bewerten wir also die Generalisierungsleistung des Modells? Eine Option besteht darin, eine Fehlermetrik auf Grundlage der Trainingsdaten zu bewerten. Leider führt dies zu verzerrten Ergebnissen, da einige Modelle auf den Trainingsdaten sehr gut abschneiden, aber nicht gut auf einem neuen Datensatz generalisieren können.
Eine zweite Methode ist die Verwendung eines Validierungsansatzes, bei dem das Trainingsset weiter in zwei Teile aufgeteilt wird: ein Trainingsset und ein Validierungsset (oder Holdout-Set). Wir können dann unser Modell auf dem neuen Trainingsset trainieren und die Leistung auf dem Validierungsset schätzen. Leider kann die Validierung mit einem einzigen Holdout-Set sehr variabel und unzuverlässig sein, es sei denn, Sie arbeiten mit sehr großen Datensätzen.
Resampling-Methoden bieten einen alternativen Ansatz, indem sie uns ermöglichen, ein Modell von Interesse wiederholt an Teilen der Trainingsdaten anzupassen und dessen Leistung auf anderen Teilen zu testen. Die beiden am häufigsten verwendeten Resampling-Methoden sind die k-fold-Kreuzvalidierung und das Bootstrap-Verfahren.
k-fold Kreuzvalidierung
Die k-Fach-Kreuzvalidierung (auch k-fold CV genannt) ist eine Resampling-Methode, bei der die Trainingsdaten zufällig in k Gruppen (auch Folds genannt) von ungefähr gleicher Größe aufgeteilt werden. Das Modell wird auf k-1 Folds angepasst und dann wird der verbleibende Fold verwendet, um die Leistung des Modells zu berechnen. Dieser Vorgang wird k Mal wiederholt; jedes Mal wird ein anderer Fold als Validierungsset behandelt. Dieser Prozess führt zu k Schätzungen des Generalisierungsfehlers.
Folglich wird bei k-Fach-CV jede Beobachtung in den Trainingsdaten einmal zurückgehalten und in das Testset einbezogen. Es gibt keine formale Regel für die Größe von k; jedoch wird bei größerem k der Unterschied zwischen der geschätzten Leistung und der tatsächlichen Leistung, die im Testset zu sehen ist, abnehmen. Andererseits kann die Verwendung eines zu großen k zu einer erhöhten Rechenlast führen.
Wir werden verschiedene Möglichkeiten behandeln, wie man Kreuzvalidierung (CV) einbeziehen kann. Wir benötigen einen Prozess, um das ML-Modell auf jede Stichprobe anzuwenden.
vfold_cv(ames, v = 10)## # 10-fold cross-validation
## # A tibble: 10 × 2
## splits id
## <list> <chr>
## 1 <split [2637/293]> Fold01
## 2 <split [2637/293]> Fold02
## 3 <split [2637/293]> Fold03
## 4 <split [2637/293]> Fold04
## 5 <split [2637/293]> Fold05
## 6 <split [2637/293]> Fold06
## 7 <split [2637/293]> Fold07
## 8 <split [2637/293]> Fold08
## 9 <split [2637/293]> Fold09
## 10 <split [2637/293]> Fold10Bootstrapping
Eine Bootstrap-Stichprobe ist eine zufällige Stichprobe der Daten, die mit Zurücklegen genommen wird. Eine Bootstrap-Stichprobe hat die gleiche Größe wie der Originaldatensatz, aus dem sie konstruiert wurde. Darüber hinaus wird die Bootstrap-Stichprobe ungefähr dieselbe Verteilung von Werten enthalten (durch Farben dargestellt) wie der Originaldatensatz.
Wir können Bootstrap-Stichproben einfach mit
rsample::bootstraps() erstellen, wie im folgenden
Codeausschnitt dargestellt.
bootstraps(ames, times = 10)## # Bootstrap sampling
## # A tibble: 10 × 2
## splits id
## <list> <chr>
## 1 <split [2930/1062]> Bootstrap01
## 2 <split [2930/1076]> Bootstrap02
## 3 <split [2930/1066]> Bootstrap03
## 4 <split [2930/1045]> Bootstrap04
## 5 <split [2930/1087]> Bootstrap05
## 6 <split [2930/1108]> Bootstrap06
## 7 <split [2930/1075]> Bootstrap07
## 8 <split [2930/1078]> Bootstrap08
## 9 <split [2930/1053]> Bootstrap09
## 10 <split [2930/1067]> Bootstrap10Bias Varianz trade-off
Vorhersagefehler können in zwei wichtige Teilkomponenten unterteilt werden: Fehler aufgrund von “Bias” und Fehler aufgrund von “Varianz”. Oft gibt es einen Kompromiss zwischen der Fähigkeit eines Modells, Bias und Varianz zu minimieren.
Bias
Bias beschreibt den Unterschied zwischen der erwarteten (oder
durchschnittlichen) Vorhersage unseres Modells und dem korrekten Wert,
den wir vorhersagen möchten. Es misst, wie weit im Allgemeinen die
Vorhersagen eines Modells vom korrekten Wert entfernt sind und gibt uns
eine Vorstellung davon, wie gut ein Modell der zugrunde liegenden
Struktur der Daten entsprechen kann.
Lineare Modelle sind klassische Beispiele für Modelle mit hoher
Verzerrung, da sie weniger flexibel sind und selten nicht-lineare,
nicht-monotone Beziehungen erfassen.
Wir müssen auch an Bias-Varianz in Bezug auf Resampling denken. Modelle mit hoher Verzerrung werden selten von dem durch das Resampling eingeführten Rauschen beeinflusst. Wenn ein Modell einen hohen Bias aufweist, wird es eine Konsistenz in seiner Resampling-Leistung aufweisen.
Varianz
Auf der anderen Seite wird der Fehler aufgrund von Varianz definiert als die Variabilität der Vorhersage eines Modells für einen bestimmten Datenpunkt. Viele Modelle (z. B. k-nearest neighbor, Entscheidungsbäume, Gradient-Boosting-Maschinen) sind sehr anpassungsfähig und bieten eine extreme Flexibilität in den Mustern, die sie anpassen können. Diese Modelle bergen jedoch ihre eigenen Probleme, da sie Gefahr laufen, sich an die Trainingsdaten anzupassen (Overfitting). Obwohl man auf den Trainingsdaten eine sehr gute Leistung erzielen kann, wird das Modell nicht automatisch gut auf unbekannte Daten generalisieren.
Da Modelle mit hoher Varianz eher zu Overfitting neigen, sind Resampling-Verfahren entscheidend, um dieses Risiko zu reduzieren.
Hyperparameter anpassen
Hyperparameter (auch Abstimmungsparameter genannt) sind die „Knöpfe“, mit denen die Komplexität von Machine-Learning-Algorithmen und somit das Bias-Varianz-Trade-off gesteuert werden können. Nicht alle Algorithmen haben Hyperparameter (z. B. die kleinsten Quadrate); die meisten haben jedoch mindestens einen oder mehrere.
Die richtige Einstellung dieser Hyperparameter hängt oft von den Daten und dem jeweiligen Problem ab und kann nicht immer allein durch die Trainingsdaten geschätzt werden. Daher benötigen wir eine Methode zur Identifizierung der optimalen Einstellung.
Als Beispiel nehmen wir das k-nearest neighbor Modell. K-nächste-Nachbarn-Modelle haben einen einzelnen Hyperparameter (k), der den vorhergesagten Wert auf Basis der k nächsten Beobachtungen im Trainingsdatensatz bestimmt. Wenn k klein ist (z. B. k = 3), macht das Modell eine Vorhersage für eine gegebene Beobachtung auf Basis des Durchschnitts der Antwortwerte für die 3 Beobachtungen im Trainingsdatensatz, die der zu vorhersagenden Beobachtung am ähnlichsten sind. Dies führt oft zu sehr variablen Vorhersagewerten, weil wir die Vorhersage (in diesem Fall einen Durchschnitt) auf eine sehr kleine Teilmenge der Trainingsdaten stützen. Wenn k größer wird, stützen wir unsere Vorhersagen auf einen Durchschnitt einer größeren Teilmenge der Trainingsdaten, was natürlich die Varianz unserer Vorhersagewerte reduziert. Kleinere k-Werte (z. B. 2, 5 oder 10) führen zu hoher Varianz (aber geringerem Bias) und größere Werte (z. B. 150) führen zu hohem Bias (aber geringerer Varianz). Der optimale k-Wert könnte irgendwo zwischen 20 und 50 liegen, aber wie wissen wir, welchen k-Wert wir verwenden sollen?
Eine Methode ist die Durchführung einer Rastersuche. Eine Rastersuche ist eine automatisierte Methode, um viele Kombinationen von Hyperparameter-Werten zu durchsuchen. Unter Verwendung von wiederholtem 10-fachem CV entspricht die Fehlerrate dem durchschnitllichen Fehler für jeden Wert von k. Für k = 46 wird der Fehler minimiert (RMSE).
Modellbewertung
Heutzutage ist es weithin akzeptiert, dass ein zuverlässigerer Ansatz zur Bewertung der Modellleistung darin besteht, die Vorhersagegenauigkeit über Verlustfunktionen zu bewerten. Verlustfunktionen sind Metriken, die die vorhergesagten Werte mit dem tatsächlichen Wert vergleichen (der Output einer Verlustfunktion wird oft als Fehler oder Pseudo-Residuum bezeichnet). Bei Resampling-Methoden bewerten wir die vorhergesagten Werte für einen Validierungssatz im Vergleich zum tatsächlichen Zielwert. Zum Beispiel kann bei Regressionen ein Weg, den Fehler zu messen, darin bestehen, die Differenz zwischen dem tatsächlichen und dem vorhergesagten Wert für eine gegebene Beobachtung zu nehmen. Der gesamte Validierungsfehler des Modells wird berechnet, indem man die Fehler über den gesamten Validierungsdatensatz aggregiert. Es gibt viele Verlustfunktionen zur Auswahl, um die Leistung eines Vorhersagemodells zu bewerten, wobei jede eine einzigartige Vorstellung von der Vorhersagegenauigkeit bietet und zwischen Regressions- und Klassifikationsmodellen variiert.
Regressionsmodelle
MSE: Mean Squared Error. \(MSE = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2.\) Ziel: minimieren
RMSE: \(RMSE = \sqrt{MSE}\). Ziel: minimieren.
Deviance: Es wird oft bei Klassifizierungsmodellen benutzt. Ziel: minimieren.
MAE: \(MAE = \frac{1}{n} \sum_{i=1}^n (|y_i - \hat{y}_i|).\) Ziel: minimieren.
RMSLE: Root mean squared logarithmic error. \(RMSLE = \sqrt{\frac{1}{n} \sum^n_{i=1}(log(y_i + 1) - log(\hat y_i + 1))^2}\). Ziel: minimieren.
\(R^2\): Ziel: maximieren.
\(R^2\) (R-Squared) und \(MSE\) (Mean Squared Error) sind beide Maße, die verwendet werden, um die Genauigkeit von Regressionsmodellen zu bewerten, aber sie messen unterschiedliche Aspekte der Vorhersageleistung. Zusammenfassend kann man sagen, dass \(MSE\) die Qualität der Vorhersagen selbst misst, während \(R^2\) die Qualität des Modells misst, um die Daten zu erklären. Ein gutes Modell hat sowohl eine niedrige \(MSE\) als auch einen hohen \(R^2\)-Wert.
Klassifizierungsmodelle
Misclassification: z.B. 13/90 falscher Klasse zugeordnet (14%). Ziel: minimieren.
Mean per class errror: Der Durchschnitt. Ziel: minimieren.
MSE: Mean squared error. Berechnet die Distanz zu 1 und quadriert. Sagt unser Modell Wahrscheinlichkeiten für A, B, C von 0.91, 0.07, 0.02 voraus, so ergibt sich für die richtige Antwort A: \(MSE = 0.09^2 = 0.0081\). Für B und C als richtige Antworten \(0.93^2\) und \(0.98^2\). Ziel: minimieren.
Log Loss/Deviance/Cross-entropy: Ziel: minimieren.
Gini Index: Ziel: minimieren.
Bei der Anwendung von Klassifikationsmodellen verwenden wir oft eine Konfusionsmatrix (confusion matrix), um bestimmte Leistungsmaße zu evaluieren. Eine confusion matrix ist einfach eine Matrix, die tatsächliche kategoriale Stufen (oder Ereignisse) mit den vorhergesagten kategorialen Stufen vergleicht.
Weitere Bewertungen können vorgenommen werden.
Accuracy: korrekt klassifiziert, z.B. \(\frac{TP+TN}{total} = \frac{100+50}{165} = 0.91\). Ziel: maximieren.
Precision: Wieviele der Vorhersagen, die wir gemacht waren, waren wirklich korrekt. \(\frac{TP}{TP+FP} = \frac{100}{100+10} = 0.91\). Ziel: maximieren.
Sensitivity: Für das eintreffende Event haben wir wieviele vorhergesagt? \(\frac{TP}{TP+FN} = \frac{100}{100+5} = 0.95\). Ziel: maximieren.
Specificity: Wie gut werden non-events klassifiziert? \(\frac{TN}{TN+FP} = \frac{50}{50+10} = 0.83\). Ziel: maximieren.
AUC: Ein guter binärer Klassifikator hat eine hohe Präzision (Precision) und Empfindlichkeit (Sensitivity). Das bedeutet, dass der Klassifikator gut abschneidet, wenn er vorhersagt, ob ein Ereignis eintritt oder nicht, was die Anzahl der falsch positiven und falsch negativen Ergebnisse minimiert. Um dieses Gleichgewicht zu erfassen, verwenden wir oft eine ROC-Kurve, die die Falsch-Positiv-Rate entlang der x-Achse und die Wahre-Positiv-Rate entlang der y-Achse darstellt. Eine Linie, die diagonal von der unteren linken Ecke zur oberen rechten Ecke verläuft, repräsentiert eine zufällige Vermutung. Je höher die Linie in der oberen linken Ecke ist, desto besser ist der Klassifikator. Der AUC berechnet den Bereich unter dieser Kurve.
Alles zusammenpacken
Geschichtetes Sampling, um unsere Daten in Trainings- und Testdatensatz zu brechen, während wir konsistente Vereteilungen zwischen den Datensätzen haben. Unser Datensatz insgesamt:
ames|>
ggplot(aes(x = Sale_Price)) +
geom_histogram(aes(y = ..density..), colour = 1) +
geom_density(lwd = 1.2, linetype = 2, colour = 2)## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Jetzt nehmen wir die Aufteilung vor. In R sind
createDataPartition() und initial_split()
Funktionen in dem Paket caret enthalten, die für die Erstellung von
Trainings- und Testdatensätzen für Machine-Learning-Modelle verwendet
werden können.
createDataPartition(): Diese Funktion teilt die Daten basierend auf einer gegebenen Prozentsatzgröße in zwei Gruppen auf (z. B. 70% für das Training und 30% für das Testen). Der Unterschied besteht jedoch darin, dass es versucht, die Klassenverteilung innerhalb der Daten zu berücksichtigen. Mit anderen Worten, es versucht, die Proportionen der Zielklasse in beiden Gruppen (Training und Test) beizubehalten.
initial_split(): Diese Funktion teilt die Daten in zwei Gruppen auf, ohne die Klassenverteilung innerhalb der Daten zu berücksichtigen. Sie können den Prozentsatz für jede Gruppe angeben, z. B. 80% für das Training und 20% für das Testen.
# Stratified sampling with the rsample package
set.seed(123)
split <- initial_split(ames, prop = 0.7,
strata = "Sale_Price")
ames_train <- training(split)
ames_test <- testing(split)Als Nächstes werden wir einen k-nearest neighbor Regressor auf unsere Daten anwenden. Dazu werden wir caret verwenden, das eine Meta-Engine ist, um die Prozesse des Resamplings, der Rastersuche und der Modellanwendung zu vereinfachen.
- Resampling: 10-fold CV, fünf mal wiederholt.
- Grid Search: Hyperparameter k = 2, 3, …, 25.
- Model training and Validation: Wir trainieren ein
knn Modell (
method = "knn") mit unserer Resampling Prozedur(trControl = cv), Grid Search (tuneGrid = hyper_grid) und unserer präferierte loss function(metric = "RMSE").
# Specify resampling strategy
cv <- trainControl(
method = "repeatedcv",
number = 10,
repeats = 5
)Hier wird ein Trainingskontroll-Objekt mit dem Namen “cv” erstellt, das als Einstellung für die Durchführung der Kreuzvalidierung beim Trainieren eines Machine-Learning-Modells verwendet werden kann.
method = “repeatedcv” gibt an, dass eine Kreuzvalidierung mit wiederholten Durchläufen (repeated cross validation) durchgeführt werden soll. Bei dieser Art der Kreuzvalidierung wird die Datensatzpartitionierung in k-Folds (hier: 10-Folds) mehrmals wiederholt, um eine robustere Schätzung der Modellgenauigkeit zu erhalten.
number = 10 gibt an, dass der Datensatz in 10 Folds aufgeteilt werden soll. In jedem Durchlauf wird ein Fold als Testdatensatz verwendet und die restlichen 9 Folds als Trainingsdatensatz verwendet.
repeats = 5 gibt an, wie oft die Kreuzvalidierung mit 10-Folds wiederholt werden soll, um eine robustere Schätzung der Modellgenauigkeit zu erhalten. Insgesamt werden also 50 Durchläufe der Kreuzvalidierung durchgeführt (10-Folds x 5 Wiederholungen).
Zusammenfassend führt das Erstellen des Trainingskontroll-Objekts mit diesen Einstellungen zu einer robusten Schätzung der Modellgenauigkeit, indem die Kreuzvalidierung mit wiederholten Durchläufen durchgeführt wird. Dadurch können die hyperparameter des Machine-Learning-Modells optimiert werden, indem verschiedene Kombinationen ausprobiert werden und diejenige mit der besten Kreuzvalidierungsgenauigkeit ausgewählt wird.
# Create grid of hyperparameter values
hyper_grid <- expand.grid(k = seq(2, 25, by = 1))# Tune a knn model using grid search
library(caret)
knn_fit <- train(
Sale_Price ~ .,
data = ames_train,
method = "knn",
trControl = cv,
tuneGrid = hyper_grid,
metric = "RMSE"
)train ist eine Funktion in der R-Package “caret” und wird verwendet, um Machine-Learning-Modelle zu trainieren. Die Funktion verwendet die übergebenen Daten, um ein Modell mit den übergebenen Methoden zu erstellen, wobei eine Kreuzvalidierung zur Modellauswahl und Hyperparameter-Optimierung verwendet wird.
Hier wird train() verwendet, um ein k-NN-Modell zu trainieren. Dazu werden folgende Parameter übergeben:
Sale_Price ~ .: Formel, die angibt, dass Sale_Price die abhängige Variable ist und alle anderen Variablen in ames_train die unabhängigen Variablen sind. data = ames_train: Der Trainingsdatensatz, der verwendet wird, um das Modell zu trainieren. method = “knn”: Gibt an, dass ein k-NN-Modell verwendet wird. trControl = cv: Das zuvor erstellte Trainingskontroll-Objekt, das angibt, wie die Kreuzvalidierung durchgeführt werden soll. tuneGrid = hyper_grid: Das zuvor erstellte Dataframe, das alle möglichen Werte für den hyperparameter k enthält. metric = “RMSE”: Gibt an, dass die Root-Mean-Square-Error (RMSE) als Bewertungsmetrik für die Kreuzvalidierung verwendet werden soll.
Die Funktion führt nun eine Kreuzvalidierung mit den angegebenen Einstellungen durch, indem sie das k-NN-Modell mit verschiedenen Werten von k trainiert und bewertet. Der beste Wert für k, der das beste Ergebnis bei der Kreuzvalidierung liefert, wird verwendet, um das endgültige Modell zu trainieren. Das trainierte Modell wird in der Variablen knn_fit gespeichert.
knn_fit
ggplot(knn_fit)Wenn wir uns die Ergebnisse ansehen, stellen wir fest, dass das beste Modell mit k = 6 übereinstimmt, was zu einem RMSE von 43.846,05 führt. Dies bedeutet, dass unser Modell im Durchschnitt den erwarteten Verkaufspreis eines Hauses um 43.85 Dollar falsch vorhersagt.
Wir haben möglicherweise das optimale k-nearest neighbor-Modell für unseren Datensatz identifiziert, aber das bedeutet nicht, dass wir das insgesamt beste mögliche Modell gefunden haben.
Feature and Target Engineering
Feature-Engineering und Target-Engineering beziehen sich auf den Prozess der Manipulation und Transformation von Merkmalen (Features) und dem Zielvariablen (Target) in einem Datensatz, um die Leistung von Vorhersagemodellen zu verbessern. Feature-Engineering beinhaltet die Identifikation, Extraktion und Umwandlung von Merkmalen, um sie für die Verwendung in Machine-Learning-Modellen besser geeignet zu machen. Target-Engineering bezieht sich auf die Modifikation des Zielmerkmals, um es besser für die Vorhersage zu machen. Beide Prozesse können dazu beitragen, die Genauigkeit und Leistung von Vorhersagemodellen zu verbessern. Die Zeit, die für die Identifizierung von Datenverarbeitungsanforderungen aufgewendet wird, kann erheblich sein und erfordert, dass Sie erhebliche Zeit damit verbringen, Ihre Daten zu verstehen.
Voraussetzungen
library(visdat)## Warning: Paket 'visdat' wurde unter R Version 4.2.3 erstelltlibrary(recipes) ## Warning: Paket 'recipes' wurde unter R Version 4.2.3 erstellt##
## Attache Paket: 'recipes'## Das folgende Objekt ist maskiert 'package:stats':
##
## stepZielvariable
Die Transformation der Zielvariable kann zu einer Verbesserung der Vorhersage führen. Schauen wir uns bei der einfachen linearen Regression die Residuen an. Diese können z.B. heavy tails haben und rechtsschief sein.
ames <- AmesHousing::make_ames()
mdl <- lm(Sale_Price~Year_Built, data = ames)
log.mdl <- lm(log1p(Sale_Price) ~ log1p(Year_Built), data = ames)ames|>
ggplot(aes(x = Year_Built, y = Sale_Price)) +
geom_point()
ames|>
ggplot(aes(x = Year_Built, y = Sale_Price)) +
geom_point() +
scale_x_log10() + scale_y_log10()
ggplot(mdl, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, color = "red")
ggplot(log.mdl, aes(x = .fitted, y = .resid)) +
geom_point() +
geom_hline(yintercept = 0, color = "red")
ames|>
ggplot(aes(x = mdl$residuals)) +
geom_histogram(fill = "steelblue", color = "black") +
labs(title = "Histogram of Residuals", x = "Residuals", y = "Frequency")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ames|>
ggplot(aes(x = log.mdl$residuals)) +
geom_histogram(fill = "steelblue", color = "black") +
labs(title = "Histogram of Residuals", x = "Residuals", y = "Frequency")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Das einfache Modell \(\text{Sale_Price}=\beta_{0} + \beta_{1} \text{Year_Built} + \epsilon\) einen nimmt normalverteilten Fehlerterm als Annahme an. Eine log Transformation kann schon hier helfen diese Annahme zu erfüllen.
Es gibt zwei Hauptansätze, um zu helfen, positiv schief verteilte Zielvariablen zu korrigieren.
Normalisierung Die meisten rechtsschiefen Verteilungen können so transformiert werden, dass sie normalverteilt sind. Transformiere dafür den Trainings- und Testdatensatz.
Benutze eine Box Cox Transformation.
Umgang mit fehlenden Werten
Features filtern
numerische Merkmalskonstruktion
Dies bezieht sich auf den Prozess der Erstellung neuer numerischer Merkmale aus vorhandenen Daten, um die Leistung von maschinellen Lernmodellen zu verbessern. Dabei können z.B. vorhandene numerische Variablen kombiniert oder transformiert werden, um neue Merkmale zu erzeugen, die einen besseren Einblick in die zugrunde liegenden Muster und Zusammenhänge der Daten bieten können.
kategoriale Merkmalskonstruktion
Dabei handelt es sich um den Prozess, neue kategoriale Merkmale aus vorhandenen Daten zu erstellen, um die Leistung von maschinellen Lernmodellen zu verbessern. Hierbei können z.B. bestehende kategoriale Variablen neu gruppiert, transformiert oder umkodiert werden, um neue Merkmale zu erzeugen, die ein besseres Verständnis der zugrunde liegenden Muster und Zusammenhänge der Daten ermöglichen.
Dimensionsreduktion
korrekte/angemessene Umsetzung
Lineare Regression
Die lineare Regression ist ein grundlegendes statistisches Modell und einer der einfachsten Algorithmen für überwachtes Lernen. Obwohl es im Vergleich zu einigen der moderneren statistischen Lernansätze, die in späteren Kapiteln beschrieben werden, etwas langweilig erscheinen mag, ist die lineare Regression immer noch eine nützliche und weit verbreitete statistische Lernmethode. Außerdem dient sie als guter Ausgangspunkt für fortschrittlichere Ansätze. Wie wir in späteren Kapiteln sehen werden, können viele der anspruchsvolleren statistischen Lernmethoden als Verallgemeinerungen oder Erweiterungen der einfachen linearen Regression betrachtet werden. Daher ist es wichtig, ein gutes Verständnis der linearen Regression zu haben, bevor man komplexere Lernmethoden studiert.
Voraussetzungen
# Modeling packages
library(caret) # for cross-validation, etc.
# Model interpretability packages
library(vip)## Warning: Paket 'vip' wurde unter R Version 4.2.3 erstellt##
## Attache Paket: 'vip'## Das folgende Objekt ist maskiert 'package:utils':
##
## viSchätzung der unbekannten Parameter
Es gibt mehrere Möglichkeiten, um “beste Anpassung” zu messen, aber das LS-Kriterium findet die “beste Anpassung” durch Minimierung der Summe der quadratischen Abweichungen (RSS):
\(RSS\left(\beta_0, \beta_1\right) = \sum_{i=1}^n\left(Y_i - \beta_0 - \beta_1 X_i\right)^2.\)
In unserem Ames housing Datensatz können wir leicht eine lineare Beziehung herstellen zwischen der Gesamtwohnfläche und dem Verkaufspreis.
model1 <- lm(Sale_Price ~ Gr_Liv_Area, data = ames_train)Unser einfaches lineares Regressionsmodell/OLS hat also einen y-Achsenabschnitt und einen Steigungsparameter.
\[\begin{equation} \widehat{Y}_{new} = \widehat{\beta}_0 + \widehat{\beta}_1 X_{new}, \end{equation}\]
options(scipen = 999)
ames_train|>
mutate(predicted = predict(model1), residuals = residuals(model1)) |>
ggplot(aes(x = Gr_Liv_Area, y = Sale_Price)) +
xlim(0, 6000) +
ylim(0, 800000) +
geom_smooth(method = 'lm', se = FALSE, color = "lightgrey") +
geom_segment(aes(xend = Gr_Liv_Area, yend = predicted), alpha = .2) +
geom_point(aes(color = abs(residuals))) +
scale_color_continuous(low = "black", high = "red") +
scale_y_continuous(labels = scales::comma) +
guides(color = FALSE) +
geom_point(aes(y = predicted), shape = 1) +
ggtitle("Regressionsmodell mit \ngeschätzter Gerade und Residuen") +
theme_bw() ## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using formula = 'y ~ x'
Um einen deteillierten Report zu erhalten, können wir
summary() benutzen.
summary(model1)##
## Call:
## lm(formula = Sale_Price ~ Gr_Liv_Area, data = ames_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -474682 -30794 -1678 23353 328183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15938.173 3851.853 4.138 0.0000365 ***
## Gr_Liv_Area 109.667 2.421 45.303 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56790 on 2047 degrees of freedom
## Multiple R-squared: 0.5007, Adjusted R-squared: 0.5004
## F-statistic: 2052 on 1 and 2047 DF, p-value: < 0.00000000000000022Der durchschnittliche Verkaufspreis steigt also pro Quadratfuß um
114.88 Einheiten.
Die Residuen sind die Differenz zwischen den beobachteten Response
Values und den vorhergesagten. Die Quantile geben einen schnellen
Überblick über ihre Verteiung. Der Standardfehler bei einer einfachen
linearen Regression misst die durchschnittliche Abweichung der
Koeffizientenschätzung für die Steigung von der tatsächlichen wahren
Steigung. In diesem Beispiel ist der Standardfehler für die Steigung des
Regressionsmodells 2.448. Das bedeutet, dass wir erwarten würden, dass
unsere Schätzung der Steigung um etwa 2.448 vom wahren Wert abweicht.
Mit anderen Worten, wenn wir das Modell mehrmals ausführen würden,
würden wir erwarten, dass die Schätzung der Steigung jedes Mal um etwa
diesen Betrag variiert. Der Koeffizient t-Wert ist ein Maß dafür, wie
viele Standardabweichungen unsere Koeffizientenschätzung von 0 entfernt
ist. Wir möchten, dass er weit von Null entfernt ist, da dies darauf
hinweist, dass wir die Nullhypothese ablehnen können.
Ein kleiner p-Wert deutet darauf hin, dass es unwahrscheinlich ist, dass wir eine Beziehung zwischen den Prädiktor und der Response Variablen aufgrund des Zufalls beobachten werden. Typischerweise ist ein p-Wert von 5% oder weniger ein guter Grenzwert.
Der Standardfehler der Residuen ist ein Maß für die Qualität der Anpassung einer linearen Regression. Theoretisch wird angenommen, dass jedes lineare Modell einen Fehlerterm E enthält. Aufgrund des Vorhandenseins dieses Fehlerterms sind wir nicht in der Lage, unsere abhängige Variable aus unserer unabhängigen Variable perfekt vorherzusagen. Der Standardfehler der Residuen ist der durchschnittliche Betrag, um den die abhängige Variable von der wahren Regressionslinie abweicht. Der Standardfehler der Residuen berechnet sich durch die Wurzel der Varianz:
\[\begin{equation} \widehat{\sigma}^2 = \frac{1}{n - p}\sum_{i = 1} ^ n r_i ^ 2, \end{equation}\]
Die Differenz zwischen den beobachteten und den vorhergesagten y-Werten wird quadriert und aufsummiert. Freiheitsgrade sind n-p, also der Stichprobenumfang weniger die Anzahl der geschätzten Parameter.
\(\widehat{\sigma}^2\) ist der MSE
und die Wurzel der RMSE, welcher in R durch sigma()
berechnet werden kann.
sigma(model1)## [1] 56787.94sigma(model1)^2## [1] 3224869786Das Konfidenzintervall für die Koeffizienten erhalten wir durch:
\[\begin{equation} \widehat{\beta}_j \pm t_{1 - \alpha / 2, n - p} \widehat{SE}\left(\widehat{\beta}_j\right). \tag{4.3} \end{equation}\]
In R durch:
confint(model1, level = 0.95)## 2.5 % 97.5 %
## (Intercept) 8384.213 23492.1336
## Gr_Liv_Area 104.920 114.4149Die R-Quadrat Statistik gibt an, wie gut das Modell die tatsächlichen Daten abbildet. Sie wird als Anteil der Varianz angegeben und misst die lineare Beziehung zwischen der unabhängigen Variable und der abhängigen Variable. Der Wert von \(R^2\) liegt immer zwischen 0 und 1. Eine Zahl nahe bei 0 zeigt an, dass die Regression die Varianz der abhängigen Variable nicht gut erklärt, während eine Zahl nahe bei 1 bedeutet, dass die beobachtete Varianz der abhängigen Variable gut erklärt wird.
Multiple lineare Regression
In der Realität haben wir natürlich oft mehr als einen Prädiktor. Hängt zum Beispiel Sale_Price ab von Gr_Liv_Area und Year_Built. Mit den Prädiktoren: Jahr, in dem das Haus gebaut wurde und Bodenfläche in Quadratfuß:
\[\begin{equation} Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon, \end{equation}\]
wobei \(X_1\) und \(X_2\) die interessierenden Features sind:
\(X_1\) Gr_Liv_Area und
\(X_2\) Year_Built.
(model2 <- lm(Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames_train))##
## Call:
## lm(formula = Sale_Price ~ Gr_Liv_Area + Year_Built, data = ames_train)
##
## Coefficients:
## (Intercept) Gr_Liv_Area Year_Built
## -2102904.57 93.83 1086.85Wir können einen Contour Plot erstellen, der uns die Haupteffekte
anzeigt. Die Oberfläche ist flach, da nur Haupteffekte vorliegen. Hängt
der Effekt eines Prediktors von dem Auftreten anderer Prediktoren ab, so
liegt ein Interaction Effect vor. Sie können durch Produkte von Features
dargesellt werden. Eine Interaktion zwischen $X_1 = $
Gr_Liv_Area und $X_2 = $ Year_Built kann also
wie folgt drgestellt werden:
\[\begin{equation} Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1 X_2 + \epsilon. \end{equation}\]
In R kann eine Interaktion durch Doppelpunkte durchgeführt werden.
lm(Sale_Price ~ Gr_Liv_Area + Year_Built + Gr_Liv_Area:Year_Built, data = ames_train)##
## Call:
## lm(formula = Sale_Price ~ Gr_Liv_Area + Year_Built + Gr_Liv_Area:Year_Built,
## data = ames_train)
##
## Coefficients:
## (Intercept) Gr_Liv_Area Year_Built
## -810498.3134 -728.5084 430.8755
## Gr_Liv_Area:Year_Built
## 0.4168Zwei Contour Plots zeigen die Unterschiede zwischen “mit” und “ohne” Interaktionseffekt auf.
Ein Modell mit nur Haupteffekten:
model3 <- lm(Sale_Price ~ ., data = ames_train)
broom::tidy(model3) ## # A tibble: 300 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.73e7 11016450. -2.47 0.0134
## 2 MS_SubClassOne_Story_1945_and_Older 3.90e3 3575. 1.09 0.275
## 3 MS_SubClassOne_Story_with_Finished_Atti… -5.39e3 12164. -0.443 0.658
## 4 MS_SubClassOne_and_Half_Story_Unfinishe… -4.41e2 13942. -0.0316 0.975
## 5 MS_SubClassOne_and_Half_Story_Finished_… 1.04e3 7250. 0.143 0.886
## 6 MS_SubClassTwo_Story_1946_and_Newer -6.67e3 5510. -1.21 0.226
## 7 MS_SubClassTwo_Story_1945_and_Older 1.57e3 6074. 0.259 0.795
## 8 MS_SubClassTwo_and_Half_Story_All_Ages 3.41e3 10149. 0.336 0.737
## 9 MS_SubClassSplit_or_Multilevel -6.67e3 11673. -0.571 0.568
## 10 MS_SubClassSplit_Foyer 1.49e3 7512. 0.199 0.843
## # ℹ 290 more rowsBewertung der Modelle
Wir haben jetzt drei verschiedene Modelle: ein Modell mit einem Prädiktor, ein Modell mit zwei Features und ein Modell, in dem alle Varianblen mit aufgenommen wurden.
Wir benutzen die RMSE Metrik und Cross-Validation, um das beste
Modell zu bestimmen. Wir können die Funktion caret::train()
verwenden, um ein lineares Modell (d.h. method="lm") mit
Kreuzvalidierung zu trainieren. Der Vorteil von caret
besteht darin, dass es eingebaute Kreuzvalidierungsfunktionen bietet,
während die Funktion lm() dies nicht tut. Der folgende
Codeabschnitt verwendet caret::train(), um Modell1 mit 10-facher
Kreuzvalidierung neu anzupassen:
set.seed(123) # for reproducibility
(cv_model1 <- train(
form = Sale_Price ~ Gr_Liv_Area,
data = ames_train,
method = "lm",
trControl = trainControl(method = "cv", number = 10)
))## Linear Regression
##
## 2049 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 1843, 1844, 1844, 1844, 1844, 1844, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 56644.76 0.510273 38851.99
##
## Tuning parameter 'intercept' was held constant at a value of TRUEWenn es auf unbekannte Daten angewendet wird, sind die Vorhersagen, die dieses Modell macht, im Durchschnitt etwa $56.410,89 vom tatsächlichen Verkaufspreis entfernt.
Schauen wir uns die anderen Modelle an:
# model 2 CV
set.seed(123)
cv_model2 <- train(
Sale_Price ~ Gr_Liv_Area + Year_Built,
data = ames_train,
method = "lm",
trControl = trainControl(method = "cv", number = 10)
)
# model 3 CV
set.seed(123)
cv_model3 <- train(
Sale_Price ~ .,
data = ames_train,
method = "lm",
trControl = trainControl(method = "cv", number = 10)
)## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschensummary(resamples(list(
model1 = cv_model1,
model2 = cv_model2,
model3 = cv_model3
)))##
## Call:
## summary.resamples(object = resamples(list(model1 = cv_model1, model2
## = cv_model2, model3 = cv_model3)))
##
## Models: model1, model2, model3
## Number of resamples: 10
##
## MAE
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## model1 34076.73 37656.23 39785.18 38851.99 40200.92 42058.68 0
## model2 29227.14 30885.17 32003.59 31695.48 32710.41 33942.26 0
## model3 14740.86 15377.88 17564.13 17559.41 19344.23 21180.02 0
##
## RMSE
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## model1 45604.65 55896.58 57000.74 56644.76 59544.08 66198.59 0
## model2 37174.26 42650.00 46869.84 46865.68 51155.14 55780.47 0
## model3 20737.09 24858.60 40515.19 41691.74 55969.21 69879.47 0
##
## Rsquared
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## model1 0.4230788 0.4621034 0.5090642 0.5102730 0.5681246 0.5996400 0
## model2 0.5829425 0.6075293 0.6865871 0.6631008 0.6976664 0.7254572 0
## model3 0.5241423 0.6343368 0.7483874 0.7578796 0.9137443 0.9328876 0Wenn wir die Ergebnisse für jedes Modell betrachten, stellen wir fest, dass wir durch das Hinzufügen weiterer Informationen über mehr Prädiktoren die Metriken zur Leistungsbewertung der Kreuzvalidierung für Out-of-Sample verbessern können. In dem Fall erzielt das Modell mit allen möglichen Haupteffekten die “beste” Leistung (im Vergleich zu den anderen beiden Modellen).
Model - Bedenken
Es gibt einige Annahmen bei der linearen Regression, die oftmals verletzt werden. Die Aufnahme von weiteren Variablen ist nicht zweifelsfrei zu empfehlen. Eine Verletzung der Annahmen kann zu fehlerhaften Interpretationen und Vorhersagen führen.
1.Lineare Beziehung Wir gehen von einer linearen Beziehung zwischen den unabhängigen und der abhängigen Variablen aus. Nicht-lineare Beziehungen können jedoch durch Transformationen der a oder/und ua Variablen in lineare Beziehungen umgewandelt werden. Durch logarithmieren zum Beispiel. Hier wird y durch log transformiert.
p1 <- ggplot(ames_train, aes(Year_Built, Sale_Price)) +
geom_point(size = 1, alpha = .4) +
geom_smooth(se = FALSE) +
scale_y_continuous("Sale price", labels = scales::dollar) +
xlab("Year built") +
ggtitle(paste("Non-transformed variables with a\n",
"non-linear relationship."))
p2 <- ggplot(ames_train, aes(Year_Built, Sale_Price)) +
geom_point(size = 1, alpha = .4) +
geom_smooth(method = "lm", se = FALSE) +
scale_y_log10("Sale price", labels = scales::dollar,
breaks = seq(0, 400000, by = 100000)) +
xlab("Year built") +
ggtitle(paste("Transforming variables can provide a\n",
"near-linear relationship."))
gridExtra::grid.arrange(p1, p2, nrow = 1)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using formula = 'y ~ x'
2.Konstante Varianz der Residuen Wenn die Fehlervarianz nicht konstant ist, sind die p-Werte und Konfidenzintervalle für die Koeffizienten ungültig. Ähnlich wie bei der Annahme einer linearen Beziehung kann eine nicht konstante Varianz oft durch Variablentransformationen oder durch Hinzufügen weiterer Prädiktoren gelöst werden. Im Folgenden sehen wir die standardisierten Residuen vs die vorhergesagten Werte unseres ersten und dritten Modells. Nur ein Modell hat konstante Varianz.
ggplot(data.frame(fitted(model1), residuals(model1)), aes(fitted(model1), residuals(model1))) + geom_point()
df1 <- broom::augment(cv_model1$finalModel, data = ames_train)
p1 <- ggplot(df1, aes(.fitted, .std.resid)) +
geom_point(size = 1, alpha = .4) +
xlab("Predicted values") +
ylab("Residuals") +
ggtitle("Model 1", subtitle = "Sale_Price ~ Gr_Liv_Area")
df2 <- broom::augment(cv_model3$finalModel, data = ames_train)
p2 <- ggplot(df2, aes(.fitted, .std.resid)) +
geom_point(size = 1, alpha = .4) +
xlab("Predicted values") +
ylab("Residuals") +
ggtitle("Model 3", subtitle = "Sale_Price ~ .")
gridExtra::grid.arrange(p1, p2, nrow = 1)## Warning: Removed 20 rows containing missing values (`geom_point()`).
3.Keine Autokorrelation Die lineare Regression geht davon aus, dass die Fehler unabhängig und unkorreliert sind. Es besteht ein deutliches Muster, das darauf hindeutet, dass Informationen über ϵ1 Informationen über ϵ2 liefern. Dieses Muster resultiert daraus, dass die Daten nach Nachbarschaft sortiert sind, was wir in diesem Modell nicht berücksichtigt haben. Folglich sind die Residuen für Häuser in derselben Nachbarschaft korreliert (Häuser in einer Nachbarschaft sind in der Regel gleich groß und können oft ähnliche Merkmale aufweisen). Da der Nachbarschaftsprädiktor in Modell 3 enthalten ist (rechtes Diagramm), wird die Korrelation in den Fehlern reduziert.
df1 <- mutate(df1, id = row_number())
df2 <- mutate(df2, id = row_number())
p1 <- ggplot(df1, aes(id, .std.resid)) +
geom_point(size = 1, alpha = .4) +
xlab("Row ID") +
ylab("Residuals") +
ggtitle("Model 1", subtitle = "Correlated residuals.")
p2 <- ggplot(df2, aes(id, .std.resid)) +
geom_point(size = 1, alpha = .4) +
xlab("Row ID") +
ylab("Residuals") +
ggtitle("Model 3", subtitle = "Uncorrelated residuals.")
gridExtra::grid.arrange(p1, p2, nrow = 1)## Warning: Removed 20 rows containing missing values (`geom_point()`).
4.Keine Multikollinearität Multikollinearität bezieht sich auf die Situation, in der zwei oder mehr Prädiktorvariablen eng miteinander verbunden sind. Das Vorhandensein von Multikollinearität kann bei der OLS-Analyse Probleme verursachen, da es schwierig sein kann, die individuellen Effekte kollinearer Variablen auf die Antwortvariable zu trennen. Tatsächlich kann Multikollinearität dazu führen, dass Prädiktorvariablen statistisch insignifikant erscheinen, obwohl sie tatsächlich signifikant sind. Dies führt offensichtlich zu einer ungenauen Interpretation der Koeffizienten und erschwert die Identifizierung einflussreicher Prädiktoren.
In Ames haben beispielsweise die Variablen Garage_Area und Garage_Cars eine Korrelation von 0,89, und beide Variablen sind eng mit unserer Antwortvariablen (Sale_Price) verbunden. Wenn wir unser vollständiges Modell betrachten, in dem beide Variablen enthalten sind, sehen wir, dass Garage_Cars statistisch signifikant ist, aber Garage_Area nicht:
summary(cv_model3) |>
broom::tidy() |>
filter(term %in% c("Garage_Area", "Garage_Cars"))## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Garage_Cars 3151. 1733. 1.82 0.0693
## 2 Garage_Area 11.9 5.69 2.10 0.0359Wenn wir ohne Garage_Cars refitten, dann ist
Garage_Area hoch signifikant.
set.seed(123)
mod_wo_Garage_Cars <- train(
Sale_Price ~ .,
data = select(ames_train, -Garage_Cars),
method = "lm",
trControl = trainControl(method = "cv", number = 10)
)## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschen
## Warning in predict.lm(modelFit, newdata): Vorhersage durch Fit ohne vollen Rang
## mag täuschensummary(mod_wo_Garage_Cars) |>
broom::tidy() |>
filter(term == "Garage_Area")## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Garage_Area 19.4 3.96 4.89 0.00000109Dies spiegelt die Instabilität im linearen Regressionsmodell wider, die durch Beziehungen zwischen Prädiktoren verursacht wird.
Darüber hinaus kann Multikollinearität auftreten, wenn ein Merkmal linear mit zwei oder mehr Merkmalen zusammenhängt (was schwieriger zu erkennen ist). In diesen Fällen ist eine manuelle Entfernung bestimmter Prädiktoren möglicherweise nicht möglich. Infolgedessen bietet der folgende Abschnitt zwei einfache Erweiterungen der linearen Regression an, bei denen die Dimensionsreduzierung vor der Durchführung der linearen Regression angewendet wird.
Principal Component Regression (PCR)
Es kann die Hauptkomponentenanalyse verwendet werden, um korrelierte Variablen mit einer geringeren Anzahl unkorrelierter Merkmale (genannt Hauptkomponenten) darzustellen, und die resultierenden Komponenten können als Prädiktoren in einem linearen Regressionsmodell verwendet werden. Dieser zweistufige Prozess wird als Hauptkomponentenregression (PCR) bezeichnet.
Die Durchführung von PCR mit Caret ist eine einfache Erweiterung
unseres vorherigen Modells. Wir geben einfach
"method = "pcr" innerhalb von train() an, um
eine PCA für alle unsere numerischen Prädiktoren durchzuführen. Du
kannst einen signifikanten Rückgang des Vorhersagefehlers im Vergleich
zu unseren vorherigen linearen Modellen erkennen, indem wir nur fünf
Hauptkomponenten verwenden, gefolgt von einem allmählichen Abfall
danach. Du wirst jedoch feststellen, dass es alle 100 Hauptkomponenten
benötigt, um einen minimalen RMSE zu erreichen.
# perform 10-fold cross validation on a PCR model tuning the
# number of principal components to use as predictors from 1-100
set.seed(123)
cv_model_pcr <- train(
Sale_Price ~ .,
data = ames_train,
method = "pcr",
trControl = trainControl(method = "cv", number = 10),
preProcess = c("zv", "center", "scale"),
tuneLength = 100
)
# model with lowest RMSE
cv_model_pcr$bestTune## ncomp
## 100 100## ncomp
## 97 97
# results for model with lowest RMSE
cv_model_pcr$results %>%
dplyr::filter(ncomp == pull(cv_model_pcr$bestTune))## ncomp RMSE Rsquared MAE RMSESD RsquaredSD MAESD
## 1 100 32522.85 0.8352098 20446.33 7017.74 0.06480652 1735.195## ncomp RMSE Rsquared MAE RMSESD RsquaredSD MAESD
## 1 97 30135.51 0.8615453 20143.42 5191.887 0.03764501 1696.534
# plot cross-validated RMSE
ggplot(cv_model_pcr)
Durch die Kontrolle der Multikollinearität mit PCR können wir im Vergleich zu den zuvor erhaltenen linearen Modellen eine signifikante Verbesserung unserer Vorhersagegenauigkeit erfahren (eine Reduktion der kreuzvalidierten RMSE auf fast $30.000).
Es ist wichtig zu beachten, dass PCR ein zweistufiger Prozess ist, bei dem der PCA-Schritt keine Aspekte der Antwort berücksichtigt, wenn er die Komponenten auswählt. Es versucht lediglich, die Variabilität im Prädiktorraum zu reduzieren. Wenn diese Variabilität mit der Variabilität der Antwort zusammenhängt, hat PCR eine gute Chance, eine Vorhersagebeziehung zu identifizieren, wie in unserem Fall. Wenn die Variabilität im Prädiktorraum jedoch nicht mit der Variabilität der Antwort zusammenhängt, kann es schwierig sein, eine Vorhersagebeziehung zu identifizieren, wenn tatsächlich eine vorhanden ist. Ein alternativer Ansatz zur Reduzierung des Einflusses der Multikollinearität ist die partielle kleinste Quadrate (PLS).
Partial Least Squares (PLS)
PLS kann als überwachte Verfahren zur Dimensionsreduzierung betrachtet werden. Ähnlich wie bei PCR konstruiert diese Technik auch eine Reihe von linearen Kombinationen der Eingangsvariablen für die Regression, aber im Gegensatz zu PCR wird die Response-Variable zur Unterstützung der Konstruktion der Hauptkomponenten verwendet. Somit können wir PLS als überwachtes Verfahren zur Dimensionsreduzierung betrachten, das neue Merkmale findet, die nicht nur den größten Teil der Informationen in den ursprünglichen Merkmalen erfassen, sondern auch mit der Response-Variable zusammenhängen.
Ähnlich wie bei PCR können wir leicht ein PLS-Modell anpassen, indem
wir das Methodenargument in train() ändern. Wie bei PCR ist
die Anzahl der Hauptkomponenten, die verwendet werden sollen, ein
Abstimmungsparameter, der durch das Modell bestimmt wird, das die beste
Vorhersagegenauigkeit erzielt (in diesem Fall den RMSE minimiert).
# perform 10-fold cross validation on a PLS model tuning the
# number of principal components to use as predictors from 1-30
set.seed(123)
cv_model_pls <- train(
Sale_Price ~ .,
data = ames_train,
method = "pls",
trControl = trainControl(method = "cv", number = 10),
preProcess = c("zv", "center", "scale"),
tuneLength = 30
)
# model with lowest RMSE
cv_model_pls$bestTune## ncomp
## 3 3## ncomp
## 20 20
# results for model with lowest RMSE
cv_model_pls$results %>%
dplyr::filter(ncomp == pull(cv_model_pls$bestTune))## ncomp RMSE Rsquared MAE RMSESD RsquaredSD MAESD
## 1 3 31077.42 0.8493617 19137.51 7878.004 0.07079317 1784.467## ncomp RMSE Rsquared MAE RMSESD RsquaredSD MAESD
## 1 20 25459.51 0.8998194 16022.68 5243.478 0.04278512 1665.61
# plot cross-validated RMSE
ggplot(cv_model_pls)
Interpretation der Feature
Die Bedeutung der Variablen zielt darauf ab, diejenigen Variablen zu identifizieren, die in unserem Modell am einflussreichsten sind. Bei linearen Regressionsmodellen wird dies in der Regel durch den absoluten Wert der t-Statistik für jeden verwendeten Modellparameter gemessen. Obwohl dies einfach ist, können die Ergebnisse schwer interpretierbar sein, wenn das Modell Wechselwirkungseffekte und komplexe Transformationen enthält. Bei einem PLS-Modell kann die Bedeutung der Variablen mithilfe der gewichteten Summen der absoluten Regressionskoeffizienten berechnet werden. Die Gewichte sind eine Funktion der Reduktion des RSS (Residual Sum of Squares) über die Anzahl der PLS-Komponenten und werden separat für jede abhängige Variable berechnet. Daher werden die Beiträge der Koeffizienten proportional zur Reduktion des RSS gewichtet.
Wir können vip::vip() verwenden, um die wichtigsten Variablen zu extrahieren und zu visualisieren. Das Bedeutungsmaß wird von 100 (am wichtigsten) auf 0 (am wenigsten wichtig) normiert.
vip(cv_model_pls, num_features = 20, method = "model")## Warning: Paket 'pls' wurde unter R Version 4.2.3 erstellt##
## Attache Paket: 'pls'## Das folgende Objekt ist maskiert 'package:caret':
##
## R2## Das folgende Objekt ist maskiert 'package:stats':
##
## loadings
Partial Dependence Plots (PDPs) helfen zu veranschaulichen wie sich eine feste lineare Änderung von x in einer festen linearen Änderung von y verhält und dabei die durchschnittliche Wirkung aller anderen Merkmale im Modell berücksichtigt wird (bei linearen Modellen entspricht die Steigung des PDPs dem entsprechenden Merkmal des OLS-Koeffizienten).
pdp::partial(cv_model_pls, "Gr_Liv_Area", grid.resolution = 20, plot = TRUE)
pdp::partial(cv_model_pls, "Total_Bsmt_SF", grid.resolution = 20, plot = TRUE)
pdp::partial(cv_model_pls, "First_Flr_SF", grid.resolution = 20, plot = TRUE)
pdp::partial(cv_model_pls, "Garage_Area", grid.resolution = 20, plot = TRUE)